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Abstract 


A procedure for the finite-difference numerical solution of the 
lifting potential flow about any number of arbitrarily shaped bodies Is 
given. The solution is baaed on a technique of automatic numerical gener- 
ation of a curvilinear coordiuate system having coordinate lines coinci- 
dent with the contours of all bodies in the field, regardless of their 
shapes and number. The effects of all numerical parameters involved are 
analyzed and appropriate values are recommended. Comparisons with analytic 
solutions for single Karman-Tref f tz airfoils and a circular cylinder pair 
show excellent agreement. The paper serve'* also to Illustrate the technique 
of application of the boundary-fitted coordinate systems to the numerical 
solution of partial differential equations. 
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I. INTRODUCTION 


Numerical incompreasible potential flow solutions for bodies of arbi- 
trary shape have generally fallen into three categories: 

(1) Integral equation methods, whereby various singular solutions of 
Laplace's equation are superposed to construct a solution satisfying the 
boundary conditions of the particular problem of interest. This type of 
approach is represented by the work in references 1-0. In these methods 
singular solutions of Laplace's equation are distributed on the body sur- 
face, and perhaps also in its interior, with the body surface represented 
by quadrala teral or triangular panels. The strengths of the singularities 
are then determined such that the superposition of the onset velocity 
field and that induced by the totality of the singularities satisfies the 
condition of vanishing normal velocity at the body surface at certain points. 

This approach has been carried to a high degree of refinement 
and is presently capable of treating the flow about multiple bodies of 
arbitrary shape. This procedure obviates calculation in the entire flow 
field and involves instead the solution of a matrix equation of order 
equal to the number of points of application of the boundary condition on 
the bodies. The primary output is the surface pressure distribution on the 
bodies ‘and the resulting aerodynamic coefficients, the velocity field can 
also be obtained, but this requires the evaluation of the velocity at each 
point in the field from a summation over all the singularities involved — a 
time consuming process. The determination of streamlines, or equivalently 
the stream function, from the velocity field is still another numerical 
problem itself. 
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(2) Finite element merhods, as represented by references 9-11. Here 
the calculation Is carried out In the entire Mow field, the field being 
divided Into finite elements. The flow solution is obtained by applying 
an Integral variational principle, or other integral relations, ovet the 
aggregate of elements, which leads to a matrix solution of order equal to 
the total number of elements In the field. The solution is thus obtained 
in the entire flow field. However, not all derivatives can be made con- 
tinuous across the boundaries between the various elements. 

(3) Conformal transformation, whereby the field is transformed to one 
of simple geometry on which the solution la known (two-dimensional flow 
only). The classic Theodorsen method [12] is one of this type. A compara- 
tive discussion of earlier applications of this and other procedures is given 
In [1]. Recently Ives [13] has extended this approach to multiple bodies. 

Finite difference solutions have been severely hindered In the pas': 
by the problem of fitting curved boundaries into the computational grid. 

The use of interpolation between grid points to represent boundary condi- 
tions on a curved boundary passing through a rectangular grid may lead to 
poor application of the boundary conditions. Since finite difference 
solutions depend on continuity of derivatives, the distribution of points 
at will in the field leads to difference expressions involving large num- 
bers of points, loss of repeat patterns over the field, and hence unreason- 
ably complex computer codes. 

However, if a curvilinear coordinate system with coordinate lines 
coincident with the field boundaries can be found, these problems vanish, 
and the finite difference approach can give very smooth solutions that do 
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not lack continuity of derivatives. The potential flow solution reported 
herein is based on Just such an approach. 

The present finite-difference potential flow solution utilizes a 
method of automatic numerical generation of a general boundary-fitted 
curvilinear coordinate system having coordinate lines coincident with all 
boundaries of a general multi-connected region containing any number of 
arbitrarily shaped bodies, whic ; hcis been reported earlier (Ref. 14-16). 

The curvilinear coordinates are generated es the solutions of an 
elliptic partial differential system with Dirlchlet boundary conditions, 
one coordinate being specified to be constant on each of tV.e boundaries, 
and a distribution of the other being specified along the boundaries. No 
restrictions are placed on the shape of the boundaries, which may even be 
time-dependent, and the method is not restricted to single bodies or, in 
principle, to tyo dimensions. Coordinate lines may be concentrated as 
desired along the boundaries. Procedures have also been developed to 
control the coordinate line spacing In the field by varying the generating 
elliptic system. 

Regardless of the number or shapes of the bodies and regardless of 
the spacing of the curvilinear coordinate lines, all numerical computations, 
both to generate the coordinate system and to subsequently solve the Lap- 
lace equation for the stream function are done on a rectangular grid with 
a square mesh, i.e., in the transformed plane. Although not necessary in 
the present application, it is also possible to cause the boundary- fit ted 
coordinate system to change in time however desired and still have all 
computation done on the fixed rectangular grid with square mesh [17]. 

This allows the curvilinear coordinate system in the physical plane to 



deform with a deforming body, free surface, or any other boundary, keeping 
a coordinate line always coincident with the boundary at all times. The 
physical coordinate system has been, in effect, eliminated from the prob> 
lem, at the expense of adding two elliptic equations to the original system. 

Since the curvilinear coordinate system has coordinate lines coin- 
cident with '^hc surface contours of all bodies present, all boundary 
conditions may be expressed at grid points, and normal derivatives on the 
bodies may be represented using only finite differences between grid points 
on coordinate lines, without need of any interpolation, even chough the 
coordinate system is not orthogonal at the boundary. With this method of 
boundary-fitted coordinate system generation, the treatment of fields with 
complex boundaries and any number of bodies is not inherently more difficult 
than problems with simple geometry. 

This use of numerically-generated boundary-fitted coordinate systems 
is not peculiar to the present application to potential flow, but rather 
is applicable to the numerical solution of any partial differential system. 
Other applications presently under consideration include viscous incompres- 
sible and compressible flow, free surface flows, turbulent flows, and 
solid mechanics problems. Some examples of such applications are given 
in Ref. 17-20. Documentation of the coordinate system generation procedure 
and computer program are given in Ref. 18 and 19. The procedure is briefly 
summarized in the next section. The present paper serves also to illustrate 
the technique of application of the boundary-fitted coordinate systems in the 
numerical solution of partial differential equations. 

II. BOUNDARY-FITTED CURVILINEAR COORDINATE SYSTEM 
A. Mathematical Constructio n 

Let it be desired to transform the two-dimensional, doubly-connected 
region, D, bounded by two closed contours of arbitrary shape into a rectang- 
ular region, D*, as shown in Fig. 1. The general transformation from 


th« physical plane (x.y) to the transformed plane (C»n) Is given by S 
C ■ Cfx,y), n ■ n(x,y). Since the baalc Idea of the transformation Is to 
generate transformation functions such chat all boundaries are coincident 
vich coordinate lines, the curvilinear coordinates (Cth) are taken as 
solutions of some suitable elliptic boundary value problem with one of 
these coordinates constant on the boundaries. The choice of a suitable 
elliptic system is restricted somewhat by the need for certain maximum 
principles as discussed in Ref. 18. The system given below allows con- 
siderable control to be exercized over the spacing of the curvilinear 
coordinate lines in the field: 

* fyy - 

”xx "yy ■ 

with Dlrichlet boundary conditions, n * constant ■ on n ■ constant 
■ r\^ on T 2 I C(x,y) e. multiple valued solution with a branch or t(x,y) 
specified (but not constant) on and V 2 » The curve on the physical 
plane transforms to the lower boundary, r*, of the transformed plane. 
Similarly, P 2 transforms to T*, etc. The right and left boundaries of the 
rectangular transformed plane, T* and F*, are coincident in the physical 

s 

plane. The curve which transforms to these boundaries connects F^ and 
F 2 and determines a branch cut for the multiple valued function C(x,y). 

Thus the functions and all derivatives are continuous across this rut. 

The inhomogeneous functions P(C»n) and Q(5,n) are sums of decaying 
exponentials that allow coordinate lines to be attracted to specified 
lines and/or points in the field or on the boundaries as discussed in more 
detail in Ref. 17-19. These functions, along with the derivative transfor- 
mation relations are given in the apoendix. 

Now since it is desired to do all numerical computation in the rec- 
tangular transformed plane, it is necessary to interchange the dependent and 
independent variables in (1). Thus using the relations from the appendix, 
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ax.- - 28x.^ • - J^(x.P(C,n) x„Q(C,n)l 

(,(, nn \ n 


(2a) 


where 


“>'« - ♦ Yv„„ 


X 2 ^ y 2 
n 'n 


J^ly^P(Cin) + y^Q(^.n) 


+ y^" 


(2b) 


6 • x.x + y.y_ 


^ ■ Vn ■ *n>'5 


with the traiisformed boundary conditions, x • fj^(C.nj^) on T*, y - g, (C.Hj^) 
on r*, x - f 2 (CiH 2 ) on P^, y • g 2 (C.H 2 ) on P^. (In the present application, 
X and y are ncndimenslonalized with respect to the airfoil chord.) 

The curvilinear coordinate system so generated has a constant n-line 
coincident with each boundary in the physical plane. The ^“lines may be 
spaced in any manner desired around the boundaries by specification of the 
C boundary conditions, or equivalently by specification of (x,y) at the 
equi-spaced C-points on the and ri 2 lines of the transformed plane. 

Control of the spacing of the n-lines may be exercised by varying the 
attraction parameters in the functions P(C,n) and Q((;,n) of Eq. (2) as 
discussed in Ref. IS and 19. 

The same procedure for boundary-fitted coordinate generation may be 
applied to regions that are more than doubly connected, i.e., have more 
than two closed boundaries or, equivalently, more than one body or hole 
within a single outer boundary. One possible transformation to the 
rectangular field for any number of bodies is illustrated in Fig. 2. Here 
the n-coordinate is required to be equal to the same constant on all the 
interior boundaries, i.e., on all bodies in the field. Let all the bodies 
be connected by arbitrary cuts and, similarly, one body bo connected to 
the outer boundary by an arbitrary cut. Since the n-coordinate is equal 
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to tht saint constant on all ths bodlss, it is, of coursa, equal to that 
constant on the cuts between the bodies also. By contrast, the ^-coordinate 
is taken constant on the cut between one body a .» the outer boundary. 

Since the locations of these cuts are not specified, the specification 
of n or ^ as constant on a cut does not overspecify the elliptic problem. 

Note that all bodies except one are split into two segments. Each 
cut appears twice on the transformed field boundary of course, the two 
segments corresponding to the two "sides" of the cut in the physical plane 
and thus being re-entrant boundaries with the functions of all derivatives 
continuous thereon. Thus, x and y are specified on the portions of the 
lower boundary of the transformed field that corresponds to the bodies, 
and also on the entire upper boundary, corresponding to the outer boundary 
in the physical field. The r mainJng portions of the lower boundary and 
the entire side boundaries are re-entrant boundaries, and thus neither 
require nor allow specification of (x,y) thereon. Other arrangements 
are also possible as discussed in detail in Ref. 18 and 19, two of which 
are used below. 

Again the elliptic Dirichlet problem (2) is solved to generate the 
boundary- fit ted coordinates (C»n). All computation, both to generate the 
coordinates and subsequently to solve the partial differential system of 
interest, are again done on the rectangular field with square mesh in the 
transformed plane. 

B. Numerical Implementation 

The transformed field for a single airfoil is illustrated in Fig. 3a. 
The physical coordinates of I points describing the body surface, (x,y), 
provide the boundary conditions along the j ■ 1 line, and those of I 
points on the physical remote boundary, usually a circle of radius ten or 



XX 


sort chord!, vupply ch« boundary condition! along tha J ■ J lina of the 

tranafonnad fiald. Since the aide boundaries of the tranafortned field are 

re-entrant, correaponding to the cut in the physical plane, we have f ■ 

^ • J 

^ and fj^j^ j" f 2 .j tor all J. Note that the value.^ of x and y are not 
apecified on these «idc boundaries. All derivatives in (2) are approxi- 
mated by second-order central difference expressions and An are both 
unity by construction, the actual values of C and n being immaterial since 

I 

cancellation occurs after substitution in the transformed equations.): 


(f ) * ~ (f 

^'rij 2 ''i+i,j 

^1-1. j^ 

(3a) 

' n^ij 2 ' i,j+l 

^i,j-l) 

(3b) 

(f ) ■ f - 2f 

^55^3 i+l.J ij 

^-l»j 

(3c) 

nn ij i,j+l ij 


(3d) 

4 ^‘i+l,j+l ^i+l,j 

-1 " ^i-l,j+l ^1-1, j-1^ 

(3e) 


The resulting set of 21 (J-1) nonlinear difference equations, two for each 
point (1,J) for 1 ■ 1, 2, — , I-l and J ■ 2, 3, — , J-1, are solved by 
accelerated Gauss-Seidel (SOR) iteration using overrelaxation. The 
iteration is considered to have converged when the maximum absolute change 
on the field between iterates is Jess than a specified value. A range of 
acceleration parameters was examined, and a value of 1.85 was nearly opti- 
mum for the airfoils considered. After convergence of the solution of (2), 
the values of the coefficients a, 6, y, J, at each point of the field 
are stored for use in the solution of the stream function equation. 

The transformed field for two airfoils is illustrated in Fig. 3b. 

The physical coordinates of body //2 at points 1 ■ 1 II, those of body 

at points 1 ■ I2 13, and finally the remaining points i ■ 14 1 on body 
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#2 art input as boundary conditions on the j ■ 1 line in the transformed 
plane. The remaining points, 1 ■ (II + 1)- — (12 - 1 ) and 1 ■ (IJ 1) — 
(14 - 1 ), on the J ■ 1 line are re-entrant points corresponding to the cut 
between the bodies in the physical plane. Therefore values at these 
points are not specified, but rather the relations f^^ 4 . i " ^14 « j 

^11 + k 0 " ^14 - k 2 ^ ^ - II - 1) hold. The rest of the 

procedure is unchanged from the case of a single airfoil, except that two 

difference equations at each of the points (1,1) for i ■ (II + 1) (12 - 

1) are added to the system, so that the to'.al number of equations is now 
2KJ-1) + 2(12 - II - 1). 



III. POTENTIAL FLOW SOLUTION 


A. Laplace Equatio n 

The r.wo-dlmensional IrroCatlonal flow about any number of bodies may 
be described by the Laplace equation for the stream function, i|i: 
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^ - 0 

yy 


(4) 


with boundary conditions 

v(x,y) ■ constant on each body (5a) 

4»(x,y) • y cos8 - x «?^n6 at infinity (5b) 

v;here 9 is the angle of attack -'f the free stream relative to the positive 
x-axls. Here the stream function is nondimensionalized relative to the 
airfoil chord and the free stream velocity. When transformed to the 
curvilinear coordinate system this equation becomes (see appendix) 

aij/.j. - + J^[Q(C.n)iJ^ + P(C,n)ij;-] - 0 (6) 

S'* nn m s 

where a, 6, Yi ^nd J are given above, and the tratisformed boundary 
conditions are, for a single body, 

'l'(Ctn) - on n “ n, (i.e., on T|) (7a) 

'j'(5in) “ y(Ctn 2 ) COS0 - x(C,H 2 ) sin9 on n - H 2 (i-e., on r|) (7b) 

The uniqueness is Implied by requiring that the solution be periodic in 
-00 < 5 < 00 , < n < H 2 * 6» Y» and J are calculated during the genera- 

tion of the coordinate system. Equation (6) is approximated using second- 
order, central differences for all derivatives, and the resulting differ- 
ence equation is solved by accelerated Gauss-Seidel (SOR) iteration on the 
rectangular transformed field. 
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The solution of (6) on Che transformed field is constructed in the 
same manner that has been previously described for the solut n of (2). 

The single equation (6) replaces the two equations (2a) and (2b), and Che 
boundary conditions are given by (7). The total number of difference 
equations thus is I(J-l) for a single airfoil and I(J-l) + (12 - II - 1) 
for two airfoils. 

3. Velocity 

The velocity components are calculated from the equations u ■ 

V ■ which in the transformed plane become, from the appendix, 

u - - x^w^)/J (8a) 

Velocities in the interior of the field may be obtained from these rela- 
tions using second-order, central difference expressions for all deriva- 
tives as given by Eq. (3). 

On the body surface, ■ 0, so that these expressions reduce to 
u ■ and v ■ Also, the unit tangent vector the the body sur- 

face is given by (see appendix) 

T » (ix^ + jy^)/*^ (9) 

Then the velocity component tangential to the surface is given by 

V e V • T ■ (ux + vy )/»^ ^ ^ (10) 

t » , s s J n 

On the surface, the C-derivatives are approximated by the second-order, 
central difference expressions of Eq. (3a), as in the interior of the 
field, at all points except those on the cut, i » 1 and i = I, where 
second order, one-sided expressicus are used. Thus 
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(fJ 


(fc) 


(f 


^3,1 ^^2,1 " ^^1,1^ 

(11a) 

1-2,1 ’ ^^1-1,1 ^^1,1^ 

(lib) 


The n-derivatives on the surface are approximated at all points by one- 
sided expressions. First, second, and third order expressions as follows 
were evaluated as discussed in the next section: 


(f ) , - f . , - f . , (first order) 

r)l»i 


(12a) 


(f ). , - -r (-f. - + 4f - - 3f , .) (second order) (12b) 


(f - T (2f, / - 9f, - + I8f, . - Ilf, ,) (third order) (12c) 


C. Kutta Condition 

The value of the boundary value of tp on the body, ip^, is determined 
by imposing the Kutta condition that the flow leave the sharp trailing 
edge of an airfoil smoothly. For a cusped trailing edge (zero included 
angle) this condition requires only that the velocity approach the same 
value at the trailing edge on the upper and lower surfaces of the airfoil. 
For a trailing edge with finite included angle it is required that the 
trailing edge be a stagnation point. It was found, however, that the 
requirement that the same limit be approached at the trailing edge on the 
upper and lower surfaces was superior numerically with both types of 
trailing edges. This limit condition was also applied by Giesing [1] 
as the Kutta condition with a finite trailing edge in the potential flow 
solution using superposition of singularities. 

In the present solution the Kutta condition thus was applied by 
requiring that the value of the velocity component tangential to the body 
surface extrapolated at the trailing edge from neighboring points on the 
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upper surface be equal co that: extrapolated from neighboring points on the 
lover surface. One, two, and three point extrapolations were evaluated, 
as well as the simple requir;iment that the velocity vanish at the trailing 
edge, These applicatiouo of the Kutta condition are as follows CHere 
superscript o refers to the trailing edge, and the other superscripts to 
successively distant neighboring points on the body surface as illustrated 
in Fig. 4. These points are, of course, equi-spaced ir the transformed 
plane. ): 



V 


( 1 -) 

t 





• 2v 


( 1 -) 

t 



(first order) (13a) 
(second order) (13b) 


3v^^‘*‘^ - 3v^^'*'^ + - v^°^ - 3v^^“^ - 3v^^"^ + v[^"^ (third order) (13c) 


D. Superposition of Solutions 

Since the system to be solved is linear in tj;, the solution for a single 
airfoil at any angle of attack may be obtained by superposing three compo- 
nent solutions; (1) a solution at 0* angle of attack with no circulation, 
(2) a solution at 90“ angle of attack with no circulation, and (3) a 
solution with circulation but zero free stream velocity as done by Giesing 
[1]. These three component solutions, written ii;^^^(C,n), i * 1, 2, 3, 
each satisfy Eq. (6), with the respective boundary conditions 




( 1 ) 

1,1 


a 0 


, i 


1 1 


'I' 


( 1 ) 

i.J 


y 


i,J’ 


i 


1 1 




( 2 ) 

i.l 


(14a) 

(14b) 


0 , i a 1 1 


(15a) 
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1 1 


(ISb) 




(3) 

1,1 


1 , 1 


1 1 


(16a) 


° ^ ^ 1 (16b) 

The complete solution with a»-bitrary circul-.tion then is 

'l»(5,n;X) - cos0 + sine - (17) 

The Kutta condition is then satisfied by choosing the coefficient \ such 
that the one of Eq. (13a) - (13c) being applied is satisfied, the 
tangential velocities being given by Eq. (10) with i|) from Eq. (17), usings 
one-sided difference expression analogous to Eq. (11) for the n-derlvative. 

Thus it is only necessary to solve the system of difference equations 
three times for a given airfoil. The solution at any angle of attack may 
then be obtained without re-solvlng the difference system. 

E. Surface Pressure and Force Coefficients 

The pressure coefficient at any point in the field may be obtained 
from the velocities via the Bernoulli equation, which in the present non- 
dimensional variables is 

Cp - 1 - lv|2 . (18) 

On the body surface this becomes, from (10), 

C = 1 - -3L ^ 2 
P j2 n 

with the derivative evaluated by one of the difference expressions of (12). 

The nondimensional force on the body is given by 

F “ - ^ C n ds 
^ P - 


( 20 ) 
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where n Is the unit outward normal to the surface , and ds is an increment 
of arc length along the surface. Since nds ■ k x dr, where r is the 
position vector of points on the surface and k is the unit normal to the 

two-dimensional plane, this becomes (see appendix) 

dr 

F - - # Cp (k X ~)dc - - # Cp (jx^ - iy^)dK 

The unit vector in the direction of the free stream is a ■ i cos9 + j sin0, 
and that normal to the free stream is k x a - J cose - i sin0, so that 

« « «» «i 

♦-he lift and drag coefficients are 


“ ? * r " ^ ^p (“Jf^cose - y^sln0)dC (21a) 

^ (y^cose -x^sin0)dC (21b) 

These integrals can be evaluated by numerical quadrature using either 
the trapezoidal rule or Simpson's rule, both of vhich were evaluated during 
the course of the investigation. For the former we have 


1 

^ f dC “ J (f J + fj j) + I f^ J (first order) 

i**2 


(22a) 


while for the latter, with I odd, 
1 


1-2 


# ' « ■ 3 <'l.l 'l,l> - 3 ^^2.1 - 3 ^^3 'i,l 


F. Multiple Airfoils 

With two airfoils, the boundary condition of Eq. (5a) becomes 


i|»(x,y) « ipj^ on the surface of body //I 
<l>(x,y) = on the surface of body //2 


(23a) 

(23b) 


With reference to Fig. 3b and the discussion in the previous section on 


the coordinate system solution, these boundary conditions become, ir the 
transformed field. 
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1,1 " h 

i - 12 13 

(24a) 

1,1 ’ '^2 

i - 1 11 and i - I4 1 

(24b) 


As in the case of the coordinate system solution, the remaining portions 
of the j ■ 1 line are re-entrant boundaries, s'» that points thereon are 
treated as field points rather than boundary points. The ^;-derivatives 
at th'S surface points, II, 12, 13, and 14, on the cuts between the bodies 
are also evaluated using the one-sided expressions of Eq. (11) in the 
calculation of the velocity on the surface. 

The Kutta condition must be applied on each body. Therefore, a 
fourth component solution is added, and the four component solutions each 
satisfy Eq. (6), with the respective boundary conditions 



- 0 1 - 1 II, 14 1 

1 I 1 

- 1 1 - 12 13 


(25a) 

(25b) 

(26a) 

(26b) 

(27a) 

(27b) 



0 


1 - 1 — r 


(27c) 




(3) 

l.J 




(4) 

1,1 


I 


1 - 1 — i:, 14 — I 


(28a) 


'll 


(4) . 
1 I 


1 - 12 13 


(28b) 


<li 


(4) 

l.J 


0 


1 - 1 1 


(2Sc) 


The complete solution with arbitrary circulation about each body is 

'l>(C,n;Xj,l2) ■ 'I'^^^(5,n)cos0 + ij;^^\c,n)sine + 

+ (29) 

The Kutta condition is then satisfied by choosing the coefficients 
and X 2 such that the one of Eq. (13a) - (13c) being applied is satisfied 
on each body. This requires only the simultaneous solution of two 
linear algebraic equations. Generalizing to N bodies, it is necessary to 
solve the difference equation system N + 2 times for a given multiple 
airfoil system. The solution at any orientation of the free stream may 
then be obtained without re-solving the difference system. 
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IV. EFFECT OF NUMERICAL PARAMETERS 


An extensive study was made to determine the effects of the various 
parameters involved on the accuracy of the numerical solution^ Numerical 
results for the lift and drag coefficients, the surface pressure distribu- 
tion, and the stream function contours for two Karman-Tref f iz airfoils were 
compared with the analytic solutions (Ref. 21) using several values for 
each of the parameters that must be chosen in the numerical solution. 

The numerical parameters involved in the solution are the following: 
Coordinat e System Parameters 

I - number of points on the airfoil (number of 5”llnes). 

J - number o‘: lines surrounding the airfoil (number of n-lines). 
r^ - radius of outer boundary. 

e - convergence criterion for iteration error norms. (Iteration is 

cs 

terminated when the maximum absolute change in x and y over the 
field between iterations becomes less than e ) 

ww 

Potential Flow Solution Parameters 

e__ - convergence criterion for iteration error norms. (See note 
FF 

with E above) 

ww 

Qg - order of extrapolation used in satisfaction of Kutta condition. 
(See Eq. 13) 

- order of approximation in calculation of surface velocity used 
in satisfaction of Kutta condition. (See Eq. 12) 

- order of approximation in calculation of surface velocity used in 
surface pressure calculation. (See Eq. 12) 

Qj - order of approximation of pressure integration used in calcula- 
tion of force coefficients. (See Eq. 22) 
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Tht comparison of the numerical results with the analytic solutions 
is presented In terms of the following quantities: 

- lift coefficient error (numerical minus analytic). 

~ drag coefficient error (analytic Is zero). 

N - maximum norm of stream function error (maximum absolute differ- 

m 

cnce between numerical and analytic values over entire field). 

- Euclidean norm of relative stream function error (root-mean- 
square over entire field of difference between numerical and 
analytic values relative to analytic value). 

Finally, results for these two airfoils at the zero-lift angle and a few 
other angles are also compared with the analytic solutions. 

The parameter comparison cases and the values of the numerical para- 
meters used therein are listed In Table 1, while the comparison of results 
Is given In Tables 2 and 3. The number of Iterations and the computer time 
required (UNIVAO 1106) are given correspondingly in Tables 4 and 5. The 
zero-lift angle results and results at other angles are given in Tables 
6 and 7. The airfoil contours are shown in Figs, 5 and 6, and the 
points on the contours and the Karman-Tref f tz parameters are given in 
Tables 8-10. Typical coordinate systems are shown in Fig. 7. Only 
a portion of the field is shown In each instance.) In addition, selected 
plots of the surface pressure distribution and the stream function contours 
In comparison with the analytic solution are presented in 8-17. 

In each figure the line is the analytic solution and the symbols are the 
numerical results. Finally, the numerical results for a Liebeck laminar 
airfoil are given in Figs. 18-21, in comparison with experimental data (Ref. 22). 
A. Effect of Point Distribution on Airfoil Contour 

As expected, the results generally i^nprove as the number of points on 
the airfoil contour is increased. There are, b vever, some exceptions, 
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•nd the point distribution is Important as well. The cases cited below 
with a number of points ending in *7* have points spaced at equal angles 
on the circle In the complex plane from which the airfoil was generated, 
wirh three additional points at half, quarter, and uighth spacing added 
above and below the trailing edge. 

The pressure distribution for Airfoil tfl with 37 points on the con- 
tour agrees fairly well with the analytic solution but with some small 
deviation near the slope break on the upper surface (Fig. 8a). With 67 points 
the agreement is excellent (Fig. 8b). An Increase to 127 points gives no real 
Improvement and, in fact, a single low point appears at the slope break. The 
streamlines also show some slight deviation from the analytic curves with 
37 points on the airfoil, particularly below the airfoil, but near perfect 
agreement is obtained with 67 points (Fig. 9). The lift and drag coefficient 
errors both decrease progressively as the number of points on the airfoil 
increases from 37 to 67 to 127 points for Airfoil //I (Table 2). However, the 
Euclidean error norm of the stream function increases slightly for 127 points 
after decreasing for 67 points. These trends are the same with Airfoil il2, 
except that the case of 127 points is not quite as good as that with 67 points 
in any respect except (Table 3;. 

The accuracy deteriorated significantly for both airfoils when the above- 
mentioned additional points near the trailing edge were removed. Thus the 
case of 61 points, having only equl-angular spaced points (in the complex 
plane), docs not give as good agreement in the pressure distribution, the 
streamlines, or the force coefficients (Table 2) as does the case with the 
additional points (67 points). In fact, the force coefficients show an 
order of magnitude improvement for Airfoil #1 with the addition of these 
points near the trailing edge. 
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A •ifflilar addition of points becwean cha aqul-angular points naar the 
laadlng adga, hovavar, doas not nacaaaarlly improva tha solution. Fig. and 
10 show that for Airfoil #1 tha prassura distribution with 43 points has a 
prasaure dip naar the laadlng adga that doas not occur with 37 points or In 
tha analytical curve. This affect becomes even more pronounced as more 
points near tha laadlng edge are added (47 points). The st;reamllnes also 
show some deviation from the analytical solution In the vicinity of the 
laadlng edge. The leading edge pressure peak Is, however, better repre- 
sented with tha additional points, as Is the small pressure dip that occurs 
at the slope break on the upper surface. The overall agreement for the 
pressure distribution is somewhat better with the addition of the six 
leading edge points to 37 point contour, and this is reflected In an 
improvement in the force coefficients (order-of-magnitude improvement in 
drag). The second addition of points (43 to 47) was deleterious in all 
respects. The addition of points near the leading edge of a contour with 
67 points (total of 73 points), however, did not improve the force coeffi- 
cients, and the pressure coefficient agreement deteriorated somewhat on the 
upper surface near the leading edge and near the slope break (Fig. 8b & c) . 

With Airfoil ifZ, the addition of the points between the equi-angular 
points near the leading edge (73 points total) deepened the low pressure 
spike near the leading edge beyond the analytical curve by a large amount, 
with a resultant adverse effect on the force coefficients (Fig. 15b & c). Although 
the lift did improve somewhat in one case (Increase from 43 to 47), it appears 
that a more regular spacing of points is appropriate on the smooth por> 
cions of the contour. If the spacing is to be decreased near the leading 
edge, the variation should be smooth and gradual, especially with leading 
edges having large curvature placed at Idrge angles of attack. More closely 
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•paced points should be added near the trailing edge, however. It should 
be noted that the equi-angular spacing in the complex plane gives an auto* 
matic concentration of points near the leading edge, but this concentra- 
tion is smooth and gradual. 

B. Effect of Number of Lines Surrounding the Body 

Both the force coefficients improved somewhat with both airfoils as 
the number of lines surrounding the airfoil (n-lines) was increased, the 
trend being more marked with the thicker profile of Airfoil //I. However, 
the pressure dip at Che slope break on the upper surface of Airfoil IH was 
deepened beyond the analytical, and more of a dip appeared at the trailing edge 
(Fig. lib & c). With Airfoil 1*2 there was little noticeable effect of the decreased 
number of lines cn Che pressure distribution or the streamlines. Thirty 
lines is clearly adequate for 1% accuracy in lift with the outer boundary 
located at 10 chords. 

C. Effect of Location of Outer Boundary 

Since the boundary conditions applied on the outer boundary are 
those appropriate at infinity, the outer boundary must be sufficiently 
distant from the body for accuracy. This effect was analyzed by changing 
both the outer boundary radius and the number of lines surrounding the 
body, so that the average mesh spacing would be vinchanged. As expected, 
the force coefficients do improve as the outer boundary recedes. However, 

10 chords was adequate for 1% accuracy in lift. A decrease in the outer 
boundary radius to 5 chords produced deterioration in both force coeffi- 
cients. The use of a more distant outer boundary definitely requires an 
Increase in the number of lines surrounding the body. 

The effect of an increase in the outer boundary radius from 10 to 20 
chords was not discernable in the pressure distribution and streamlines 
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for Airfoil 92 . The effect* with Airfoil #1 (Fig. 11) were similar to those 
of an increase in the number of lines surrounding the body. This results from 
the natural tendency of the coordinate system to expand outward from the body, 
so that a simultaneous increase in the outer boundary radius and the number 
of lines surrounding the body results in a closer line spacing near the body, 
even though the average spacing is unchanged. 

D . Effect of Order of Extrapolation Used in Kutta Condition 

The use of two or one point extrapolation of the velocity to the trailing 
edge in application of the Kutta condition resulted in a progressively deeper 
low pressure spike at the trailing edge with Airfoil 91 (cf Fig. 8c & 13b). 

This spike is removed with three~polnt extrapolation. This effect was less 
pronounced with Airfoil 92 since the trailing edge was less sharp. The use 
of the requirement that the trailing edge be a stagnation point gave about 
as good results as the use of the extrapolation to che trailing edge with the 
less sharp trailing edge of Airfoil 92 (cf Fig. 8c & 13a). Three-point extra- 
polation is the most appropriate condition, as is indicated by the values 
for Airfoil 91 in Table 2. 

E . Effect of Order of Velocity Diff erence Expression Used in Kutta Condition 

A first-order expression for the velocity is not sufficiently accurate 
and produced a low pressure spike at the trailing edge with Airfoil 91 . This 
spike was removed by the use of a second-order expression (Fig. 8c). Further 
increase in order gave no improvement (Fig. 13d). Again these effects are 
less evident with the less sharp trailing edge of Airfoil 92 . 

F. Effect of Pro .r of Ve l ocity Difference Expression Used in Pressure Calculation 
Calculation 

Here again a first-order expression does not give an accurate pressure 
distribution with Airfoil 91 (Fig. lAa) . A second-order expression, however. 
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Is quite accurate, (Fig. 8c) and little further improvement occured with an 
increase to third>order (Fig. 14b). 

G. Effect of Order of Pressure Integration 

Perhaps surprisingly, Sirnson's rule (second-order) did not give as 
accurate a pressure integration as did the first-order trapezoidal rule. 
This probably is a result of the fact that although the points on the 
contour are equally spaced in the transformed plane, so that Simpson's 
rule can be applied, they are not equally spaced in the physical plane. 
Therefore, the truncation error term of the Simpson integration will be 
dependent on the rate-of-change of the physical coordinate spacing, so 
that, although the order is second, the coefficient may be high with a 
resultant loss in overall accuracy in regions of rapidly varying spacing. 

H. Recommended Values 

In view of the present results, 1% accuracy in lift can be achieved 
with the following values for the numerical parameters: 

(a) 37 points on the airfoil contour, spaced with smooth and gradual 
concentration near the leading edge and with a few more closely 
spaced points added near the trailing edge. 

(b) 30 lines surrounding the airfoil. 

(c) Outer boundary at !0 chords. 

(d) Coordinate system convergence criteria of 10“*^. 

(e) Stream function convergence criterii.. of 10“^*. 

(f) Three-point extrapolation to trailing edge in satisfaction of 
Kutta condition. 

(g) Second-order difference expression for velocity. 

(h) Trapezoidal pressure integration. 


These values are the minimum considered adequate for 1% accuracy in 
lift. The changes in force coefficients encountered with more stringent 
values are summarized in Table H for Airfoil //I. 

Items (f) - (h) are adequate for even higher accuracy. However, items 
(a) ■* (c) must be increased, and items (d) - (e) decreased for higher 
accuracy. A change of the convergence criteria to 10“^ is relatively 
inexpensive in computer time. However, computer time varies approximately 
linearly with the number of points on the airfoil and quadratically with 
the number of lines surrounding the airfoil. More points on the contour 
will obviously be required with more Irregularly shaped bodies. The 
computer time is not significantly affected by the shape with the same 
number of points. 


V. MULTIPLE- BODY SOLUTIONS 


In the present study, the numerical results for potential flow about 
two circular cylinders, aligned such that a line connecting the circle 
centers is normal to the uniform free stream velocity, are compared with 
the analytic solution for the pressure distribution on the cylinder 
surfaces (Ref. 23). The cylinder axes were separated by three radii in 
all cases. The outer boundary was located at 20 cylinder diameters from 
the mid-point between the cylinders in all cases, and there were 40 n-lines 
surrounding the bodies, except as noted. 

As noted above, there are a number of different possible arrangements 
in which the boundaries in a multiple-body field may be distributed around 
the rectangular boundary of the transformed field. Many of these arrange- 
ments are illustrated in Ref. 18-19. Three such arrangements were considered 
in the present study (Fig. 22). 

If the generating partial differential system for the curvilinear 
coordinates is simply the pair of Laplace equations (P » Q = 0 in (2)), then 
the coordinate lines have a tendency to be attracted, as it were, to convex 
portions of the boundary and repelled from concave portions. This presents 
a problem with the two-body segment arrangement shown in Fig. 22a, for, 
with the n-coordinate having the same value on both .bodies , the cut 
connecting the bodies must also be a line of the same ri“value. Conse- 
quently, a concave region develops at the intersections of this cut with 
the bodies. The resulting coordinate system is shown in Fig. 23a, and the 
results for the surface pressure distribution in Fig. 23b are grossly In 
error because of the large truncation error that occurs with the widely 
spaced coordinate lines in the regions of the cut intersections. 
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However, coordinate system control, as discussed above, can be used 
to attract the coordinate lines Into these concave regions. As noted above. 
It is possible to attract the coordinate lines to specified lines and/or 
points on the boundary and In the field. Various different combinations 
of line and point attractions were considered and representative results 
are shown in Figs. 2A and 25. 

In Fig-- 2Aa-e, the coordinate systems and pressure distributions 
are shown for five amplitudes of point attraction of the n-lines (sur- 
rounding the bodies) to the cut intersections. The attraction decay factor 
(see appendix ) was 0.5 in all cases given. Too large a decay factor is 
ineffective, while too small a factor extends the attraction beyond the 
local region of interest. The pressure distribution changes drastically 
from that from that obtained without attraction (Fig. 23) as the attraction 
amplitude is initially increased. The changes become progressively smaller 
as the amplitude increases farther, with very little change from an 
amplitude of 500 to one of 2000. (Although the best fit to the analytic 
solution occurs with an amplitude of 250, this is only a fortuitous 
circumstance, since the results are still changing with attraction amplitude 
at that point.) This progressively diminishing effect of increasing 
attraction amplitude is also evident in the force coefficients normal to the 
free stream given in Table 12 . 

Comparison of Figs. 2Ad and 24f shows the effect of decreasing the 
number of C“lines passing between the bodies. The solution is not greatly 
affected by the decrease in the number of lines except in the region of 
the cut intersections, where the angle between the (;-line emanating from 
the cut intersection and the n-line between the bodies is smaller. Too 
small an angle between coordinate lines increases the local truncation 
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error, so that it is necessary to avoid having the C-line emanating from 
the cut intersection deviate greatly from 45*. This means that the point 
spacing on the cut between the bodies should be nearly the same as that on 
the body contours in the neighborhood of the cut intersections. Table 12 
shows that an addition of four lines between the bodies, by contrast, 
had little effect.' 

This need for relatively equal point spacing on the body contours and 
the cut between is further illustrated in Figs. 25b and 25d, where, with 
fewer points on the bodies, the smaller number of lines between the bodies 
gives better results. The force coefficients results in Table 12 show a 
further deterioration as still more lines are added between the bodies. 
This figure also shows that, with this smaller number of points on the 
bodies, the increasing attraction produces little change from an ampli- 
tude of 500 to one of 1000, but a larger change then occurs at a larger 
amplitude. Comparison of the force coefficients in Table 12 also shows 
that 31 points on each body is not enough for accuracy. 

That the addition of more points on the body contour in the neigh- 
borhood of the cut intersection is inferior to uniformly spaced points is 
illustrated in Figure 26a. Here the additional closely spaced points on 
the contour have caused the ^-line angles on the contour near the cut 
intersections to become too small. Some control over these angles can be 
exercised by attracting the C“lines, as well as the n~line to the cut 
intersections, the results of which are illustrated in Figure 26b. How- 
ever, the curvature of the C~lines in the field near the cut is now quite 
large, and truncation error results therefrom. Other combinations of 
attractions were also considered, and some results therefor are given in 
Table ’-2. 


Table 12 alao gives the force coefficients for some cases with closer 
and more distant outer boundaries and cases with fewer n-lines surrounding 
the bodies. These results show only small differences with 61 points on 
the bodies. 

TVo other coordinate system configurations are shovm in Figs. 27 and 
28. That of Fig. 28 is similar to a bi-polar coordinate system and pro- 
duces a near-perfect comparison with the analytic results on one cylinder 
with 61 points on each body. With 31 points, however, the accuracy is not 
quite as good. Since this configuration has different coordinate values 
on the two bodies, cut Intersections of the type discussed above do not 
occur, and the coordinate system is much more regular near the bodies. 

However, in contrast to the previous configuration, the results on the 
other cylinder are not quite identical. This appeared to be a bi-stable 
situation, with excellent comparison occuring on the other cylinder in 
some cases, a result, perhaps, of the SOR iteration sweeping the field in 
one direction. The force coefficient of the previously considered config- 
uration with 61 points on the contour (0.6534) is within 3% of the essentially 
exact value (0.6331) obtained from Fig. 28b. The results of the configuration of 
Fig. 27 could be improved by coordinate system control as used above. 

Finally, Figs. 29 and 30 show the coordinate system and potential 
solution for a multiple airfoil consisting of two Karman-Tref f tz airfoils, 
one being positioned as a separated flap. Coordinate system control has 
been used to attract the coordinate lines to the airfoils and to the cut 


between. 
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VI. CONCLUSION 

The use of body-fitted curvilinear coordinate systems allows numerical 
potential flow solutions for fields containing any number of bodies of 
arbitrary shape to be produced by finite-difference methods essentially as 
easily as that about single simple bodies. The computer code is not 
dependent on either the number or the shapes of the bodies, so that 
different bodies can be treated by simple changes in the input. Further 
investigation of the control of r.he coordinate system is presently in 
progress with the purpose of improving the coordinate configuration in 
concave regions created by cut intersections with members of multiple- 
body combinations. Multiple-body viscous solutions are also under 
development . 



APPENDIX 


DERIVATIVES AND VECTORS IN THE TRANSFORMED PLANE 


Derivative Transformations 




f ^ “ (x_f - X fj.)/J 

y '3y'x,t ' C n n C 


' ^3t'x,y " 


- i(x f - X £ )(|^) 


V2f . 


+ - (af_- - 2gf^ + Yf )/ j2 + Qf + pf 

3^2 3y2 nn n C 


(A.l) 
(A. 2) 

(A. 3) 
(A. 4) 


Unit Vectors in the Transformed Plane 


Let n be the unit vector normal to a line of constant n and t be 

-n ,n 

the unit vector tangent to an n-line. Utilizing similar definitions for unit 
normals and tangents to ^“constant lines, there results 
Vn 


"n ■ I v^r 

VC 

| 7 £| ■ 
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(A. 7) 


T- = n. X k — (x i + y ^)//a (A. 8 ) 

-5 - n. nt 

where 1, j» and k are the unit vectors parallel to the conventional cartesian 
coordinates x, y, and z, respectively. Vector components along the tangents 
and normals to lines of constant ^ and n are obtained by dotting the desired 
vector with one of the above. For example, if F » Fj^i + F 2 j, then the 
component tangent to an n-line, F , is given by 

In 

F - F • T - + yff-)/«/Y (A. 9) 

Similarly, directional derivatives of a scalar function f in the above 
directions can be obtained from the inner product of the gradient of f, 7f, 
and the appropriate unit vector. For example, the directional derivative 
normal to an n-line is 

- n • Vf - (yf - 8f,)/J»^ (A. 10) 

3n^ .n - n C 

n 

Integral Transformation 


Let S be the closed cylindrical surface of unit depth whose perimeter 
is specified by the contour in the physical plane (Flg< Dt and whose 
outward unit normal at any point is n(x,y). (Note that this is the unit normal 
to the n-line coincident with Fj^.). Then, 



(A. 11) 


where n, is the value of n on F, , C ^ and C 

1 1 min max 

values of Z, and and y^ are evaluated along n 


are the minimum and maximum 
1 * 
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Coordinate Syste m Control 

The effecL changing the functions P(C.n) and Q(C,n) on the coordi- 
nate system is discussed in some detail in Ref. 18 and 19. One particularly 
effective procedure is to choose P and Q as exponential tcn>.s, so that the 
coordinates are generated as the solutions of 




n 

a^ sgn(^ - C^)exp(-c^|C - ^,^1) 


m 


- I 

j-1 


sgn(4 - )exp (-d^/(5i - ' n^)^) 


= P(5,n) 


XX 


yy 


n 

- E a sgn(n - n,)exp(-cjn 

i-1 ^ 


(A. 12) 
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- E 
j-1 


sgn(ri - rij)exp(-d^/(4 - + (n - 


i Q(C,n) (A. 13) 

where the positive amplitudes and decay factors are not necessarily the same 
in the two equations. Here the first terms have the effect of attracting C~lines to 
the C - lines in Eq. (A. 12), and attracting n-lines to the r, = lines in 
(A. 13). The second terms cause ^-lines to be attracted to the points (^j.Hj) 
in (A. 12), and n-lines to be attracted to the points 1" (A. 13). 

Several examples of the use of coordinate system control are given in Ref. 


18 and 19. 
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Paraneter Cooparlson Cases - Airfoils #1 and #2 
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order of approximation in calculation of surface velocity used in 
calculation (see Eq. 12). 
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TABLE 2 


Comparison wlch Analytic Solution - Airfoil Ifl 

Case ACj^ ACp N 

Points on Body 


37 

16 

-0.0087 

0.0102 

3.35-01 

2.06-03 

43 

2 

0.0059 

-0.0008 

3.29-01 

3.98-03 

47 

3 

0.0120 

-0.0212 

3.29-01 

4.31-03 

61 

6 

-0.0577 

-0.0325 

3.28-01 

4.44-03 

67 

20 

-0.0019 

0.0040 

3.37-01 

1.31-03 

73 

5 

0.0019 

-0.0046 

3.29-01 

2.07-03 

127 

1 

0.0006 

0.0017 

3.40-01 

9.96-03 

Points to ®, 

at 10 




30 

32 

-0.0080 

0.0038 

3.37-01 

1.21-03 

60 

7 

0.0032 

-0.0026 

3.78-01 

1.38-03 

Points to ®, <*> 

at 20 




30 

12 

-0.0169 

0.0078 

3.65-01 

5.65-03 

60 

8 

-0.0057 

0.0010 

3.64-01 

7.89-04 

Location 

of ®, 

Step Size 10/60 avg. 



5 

11 

0.0296 

-0.0067 

3.08-01 

2.10-03 

10 

7 

0.0032 

-0.0026 

3.78-01 

1.38-03 

Location 

of ®, 

Step Size 20/60 

avg. 



10 

32 

-0.0080 

0.0038 

3.37-01 

1. 21-03 

20 

8 

-0.0057 

0.0010 

3.64-Oi 

7.89-04 

Location 

of ®, 

Step Size 40/60 avg. 

. 


20 

12 

-0.0169 

0.0078 

3.65-01 

5.65-03 

40 

9 

-0.0100 

0.0028 

3.86-01 

4.85-04 

Location 

of ®, 

30 Points to ® 




5 

11 

0.0300 

-0.0067 

3.08-01 

2.10-03 

10 

32 

-0.0080 

0.0038 

3.37-01 

1.21-03 

20 

12 

-0.0169 

0.0078 

3.65-01 

5.65-03 


TABLE 2 (CONTINUED) 



Case 

^^L 


N. 

«2 

Location 

of ®, 60 

Points to * 




10 

7 

0.0032 

-0.0026 

3.78-01 

1.38-03 

20 

8 

-0.0057 

0.0010 

3.64-01 

7.89-04 

40 

9 

-0.0100 

0.0028 

3.86-01 

4.85-04 

80 

10 

-0.0120 

0.0041 

3.56-01 

1.11-03 

Coordinate System 

Convergence 

Criteria, 37 

Points on Body 


10-2 

13 

-0.1691 

0.0049 

3.42-01 

3.80-03 

10-3 

14 

-0.0118 

0.0108 

3.36-01 

3.14-03 

10-4 

15 

-0.0108 

0.0103 

3.35-01 

2.09-03 

10-5 

16 

-O.0087 

0.0103 

3.35-01 

2.06-03 

Coordinate System 

Convergence 

Criteria, 67 

Points on Body 


10-3 

18 

-0.0345 

0.0035 

3.38-01 

1.28-03 

10-4 

19 

-0.0051 

0.0040 

3.37-01 

1.24-03 

10-5 

20 

-0.0019 

0.0040 

3.37-01 

1.31-03 

Stream Function Convergence Criteria 



10-2 

21 

-0.0687 

-0.1066 

3.32-01 

2.14-02 

10-3 

22 . 

0.0548 

-0.0084 

3.37-01 

1.22-03 

10-4 

23 

-0.0096 

0.0047 

3.37-01 

1.26-03 

10-5 

20 

-0.0019 

0.0040 

3.37-01 

1.31-03 

Order of 

Extrapolation for Kutta Condition 



A 

24 

-0.0044 

-0.0046 

3.29-01 

1.16-02 

2 

25 

0.0027 

-0.0046 

3.29-01 

1.93-03 

3 

5 

0.0019 

-0.0046 

3.29-01 

2.07-03 

Stagnation 





Point 

26 

0.0581 

-0.0050 

3.30-01 

1.18-02 

Order of 

Velocity 

for Kutta Condition 

• 


1 

27 

-0.0814 

-0.0041 

3.28-01 

6.61-03 

2 

5 

0.0019 

-0.0046 

3.29-01 

2.07-03 

3 

28 

0.0083 

-0.0047 

3.29-01 

3.08-03 

Order of 

Velocity 

for Pressure 

Calculation 



1 

29 

0.0722 

0.0109 

3.29-01 

2.07-03 

2 

5 

0.0019 

-0.0046 

3.29-01 

2.07-03 

3 

30 

0.0084 

-0.0015 

3.29-01 

2.07-03 

Order of 

Pressure 

Integration 




1 

5 

0.0019 

-0.0046 

3.29-01 

2.07-03 

2 

31 

0.0040 

-0.0087 

3.29-01 

2.07-03 



TABLE 2 (CONCLUDED) 


Note: Case numbers correspond to those of Table 1. 

Legend : 

The comparison of the numerical results with the analytic 
solutions is presented in terms of the following quantities: 

ACj^ - lift coefficient error (numerical minus analytic) 

ACp - drag coefficient error (analytic is zero) 

N 

'■' - maximum norm of stream function error (maximum absolute 

difference between numerical and analytic values over 
entire field) 

N 2 - Euclidean norm of relative stream function error (root- 
roean-square over entire, field of difference between 
numerical and analytic values relative to analytic 
value) 



TABLE 3 


Comparison with Analytic Solution - Airfoil //2 



Case 


AC, 


N 

No 




L 

D 

00 

2 

Points on 

i Body 






43 

2 


-0.0633 

0.2135 

1.11-01 

2.96-02 

47 

3 


0.0023 

0.2487 

1.10-01 

2.31-02 

53 

4 


0.0086 

0.2580 

1.09-01 

7.21-03 

61 

6 


-0.0741 

-0.0211 

1.57-01 

6.68-03 

67 

20 


0.0475 

-0.0250 

3.13-02 

1.04-03 

73 

5 


0.0890 

0.0990 

1.11-01 

2.05-02 

127 

1 


0.0973 

-0.0197 

1.11-01 

1.71-02 

Points to 

1 OU ^ 00 

at 

10 




30 

32 


0.0433 

-0.0253 

3.93-02 

2.63-03 

60 

7 


0.0307 

-0.0217 

4.11-02 

2.02-03 

Points to 

1 00^ OD 

at 

20 




30 

12 


0.0412 

-0.0280 

4.67-02 

1.48-03 

60 

8 


0.0171 

-0.0159 

4.47-02 

2.25-03 

Location 

of *, 

Step Size 10/60 

Avg. 



5 

11 


0.0929 

-0.0423 

4 . 68-02 

6.86-03 

10 

7 


0.0307 

-0.0217 

4.11-02 

2.02-03 

Location 

of 

Step Size 20/60 Avg. 



10 

32 


0.0433 

-0.0253 

3.93-02 

2.63-03 

20 

8 


0.0171 

-0.0159 

4.47-02 

2.2.5-03 

Location 

of 

Step Size 40/60 Avg. 



20. 

12 


0.0412 

-0.0280 

4.67-02 

1.48-03 

40 

9 


0.0156 

-0.0146 

5,45-02 

2.23-03 

Location 

of ®, 

30 

Points to ® 




5 

11 


0.0929 

-0.0423 

4.68-02 

6.86-03 

10 

32 


0.0433 

-0.0253 

3.93-02 

2.63-03 

20 

12 


0.0412 

-0.0280 

4.67-02 

1.48-03 

Location 

of ®, 

60 

Points to « 




10 

7 


0.0307 

-0.0217 

4.11-02 

2.02-03 

20 

8 


0.0171 

-0.0159 

4.47-02 

2.25-03 

40 

9 


0.0156 

-0.0146 

5.45-02 

2.23-03 

80 

10 


0.0176 

-0.0151 

6.86-02 

5.22-03 


TABLE 3 CONTINUED) 



Case 

AC, 


N 

No 



L 

D 

00 

2 

Coordinate System 

Convergence 

Criteria 


• 

10-2 

17 

-0.4180 

-0.0204 

4.81-02 

6.40-03 

10-3 

18 

-0.0619 

-0.0245 

3.09-02 

1.37-03 

10-4 

19 

0.0430 

-0.0249 

3.13-02 

8.89-04 

10-5 

20 

0.0475 

-0.0250 

3.13-02 

1.04-03 

Stream Function Convergence Criteria 



10-2 

21 

0.6407 

-0.2199 

4,01-01 

3.62-02 

10-3 

22 

0.1192 

-0.0620 

9.32-02 

4.18-03 

10-4 

23 

0.0487 

-0.0255 

3.92-02 

2.20-03 

10-5 

20 

0.0475 

-0.0250 

3.13-02 

1.04-03 

Order of 

Extrapolation for Kutta Condition 



1 

24 

0.0857 

0.0989 

1.11-01 

2.05-02 

2 

25 

0.0902 

0.0991 

1.11-01 

2.05-02 

3 

5 

0.0890 

0.0990 

1.11-01 

2,05-02 

Stagnation 





Point 

26 

0.0881 

0.0990 

1.11-01 

2.05-02 

Order of 

Velocity 

for Kutta Condition 



1 

27 

0.1150 

0.0997 

1.18-01 

3.62-03 

2 

5 

0.0890 

0,0990 

1.11-01 

2.05-02 

3 

28 

0.1057 

0.0995 

1.15-01 

3.11-03 

Order of 

Velocity 

for Pressure 

Calculation 



1 

29 

0.02401 

0.1397 

1.11-01 

2.05-02 

2 

5 

0.0890 

0.0990 

1.11-01 

2.05-02 

3 

30 

0.0744 

0.0917 

1.11-01 

2.05-02 

Order of 

Pressure 

Integration 




1 

5 

0.0890 

0,0990 

1.11-01 

2.05-02 

2 

31 

-0.2298 

0.0879 

1.11-01 

2.05-02 

Note: 

Case numbers correspond to those of 

Table 1. 


Legend : 







The comparison of the numerical results with the analytic 
solutions is presented in terms of the following quantities 


TABLE 3 (CONCLUDED) 


-> lift coefficient error (numerical minus analytic) 

ACjj - drag coefficient error (analytic is zero) 

- maximum norm of stream function (maximum absolute difference 
between numerical and analytic values over entire field) 

N 2 “ Euclidean norm of relative stream function error (root-mean- 
square over entire field of difference between numerical 
and analytic values relative to analytic value ) 


TABLE 4 


Iterations and Computer Time - Airfoil 01 


Points on 

Case 

> Body 

'"cs 

'^0 

37 

16 

95 

60 

43 

2 

105 

62 

47 

3 

106 

65 

61 

6 

132 

79 

67 

20 

151 

83 

73 

5 

172 

96 

127 

1 

401 

197 

Points at 

OO ^ 00 

at 10 


30 

32 

116 

53 

60 

7 

227 

73 

Points to 

00 ^ 00 

at 20 


30 

12 

165 

58 

60 

8 

290 

74 

Location 

of 

Step Sirs 

10/60 avg 

5 

11 

123 

47 

10 

7 

227 

73 

Location 

of <*>, 

Step Size 

• 

20/60 avg 

10 

32 

116 

53 

20 

8 

290 

74 

Location 

of ®, 

Step Size 

40/60 avg 

20 

12 

165 

58 

40 

9 

360 

75 

Location 

of 

30 Points 

to “> 

5 

11 

123 

47 

10 

32 

116 

53 

20 

12 

165 

58 


^^90 ^^CIR ^"^PSI 


77 

78 

0:57 

0:56 

81 

84 

1:14 

1:09 

83 

92 

1:19 

1:11 

95 

124 

2:04 

1:52 

114 

121 

2:28 

1:59 

100 

156 

3:10 

2:35 

195 

348 

11:57 

7:46 


73 

101 

1:58 

1:39 

79 

188 

7:21 

4:06 


80 

125 

2:44 

1:49 

82 

208 

10:13 

5:44 


69 

79 

2:06 

1:27 

79 

188 

7:21 

4:06 

73 

101 

1:58 

1:39 

82 

208 

10:13 

5:44 

80 

125 

2:44 

1:49 

86 

232 

11:34 

5:54 

69 

79 

2:06 

1:27 

73 

101 

1:58 

1:39 

80 

125 

2:44 

1:49 


TABLE 4 (CONCLUDED) 



Case 


'^0 

IT 

90 


^^CS 

CTpsi 

Location of 60 

Points to * 






10 

7 

227 

73 

79 

188 

7:21 

4:06 

20 

8 

290 

74 

82 

208 

10:31 

5:44 

40 

9 

360 

75 

86 

232 

11:34 

5:54 

80 

10 

446 

78 

91 

258 

14:13 

5:01 

Coordinate 

System 

Convergence 

Criteria, 37 

Points on 

Body 


10-2 

13 

38 

59 

78 

78 

0:30 

0:57 

10-3 

14 

54 

60 

77 

78 

0:39 

0:57 

10-4 

15 

75 

60 

77 

78 

0:49 

1:16 

10-5 

16 

95 

60 

77 

78 

0:57 

0:56 

Coordinate 

System 

Convergence 

Criteria, 67 

Points on Body 


10-3 

18 


83 

114 

121 


1:59 

10-4 

19 

116 

83 

114 

121 

1:58 

1:59 

10-5 

20 

151 

83 

114 

121 

2:28 

1:59 

Stream Function Convergence Criteria 





10-2 

21 

151 

6 

25 

30 

2:28 

0:40 

10-3 

22 

151 

21 

54 

60 

2:28 

1:02 

10-4 

23 

151 

51 

87 

90 

2:28 

1:30 

10-5 

20 

151 

83 

114 

121 

2:28 

1:59 


Note //I: Case numbers correspond to those of Table 1. 

Note #2: CPU time is given in MINUTES : SECONDS . These times are 

subject to some variation depending on the operating 
conditions at the time of the run. 


Legend : 


”cs - 

I terations 

for 

”o - 

Iterations 

for 

IT 

^^90 

Iterations 

for 

IT 

CIR 

Iterations 

for 


coordinate system 

0® potential flow solution 
90® potential flow solution 
circulation potential flow solution 

flow solution 


CT^g - CPU time for coordinate system 

CTpgj - CPU time for complete potential 


TABLE S 


Iterations and Computer 

Time - 

Airfoil ft2 



Case 

IT 

CS 

'^0 

'^90 

IT 

CIR 

^*cs 

"^PSI 

Points on Body 







43 i- 

95 

53 

72 

82 

1:27 

1:10 

yi 3 

104 

57 

77 

92 

1:17 

1:09 

53 ^ 

118 

59 

81 

106 

1:38 

1:21 

61 6 

144 

65 

88 

124 

2:84 

2:09 

67 20 

162 

68 

92 

141 

2:40 

2:11 

73 5 

182 

75 

100 

157 

3:25 

2:28 

127 1 

419 

131 

170 

350 

12:56 

1:01 

Points to ®, ® 

at 10 






30 32 

125 

36 

61 

101 

2:11 

2:01 

60 7 

239 

42 

75 

188 

7:54 

3:45 

Points to ®, ® 

at 20 






30 12 

225 

37 

64 

125 

4:01 

2:04 

60 8 

299 

39 

71 

208 

9:58 

3:54 

Location of ®, 

Step Size 

10/60 Avg. 





5 11 

86 

34 

57 

79 

1:33 

1:19 

10 7 

239 

42 

75 

188 

7:54 

3:45 

Location of ®, 

Step Size 

20/60 Avg. 





10 32 

125 

36 

61 

101 

2:11 

2:01 

20 8 

299 

39 

71 

208 

9:58 • 

3:54 

Location of ®, 

Step Size 

40/60 Avg. 





20 • 12 

225 

37 

64 

125 

4:01 

2:04 

40 9 

374 

37 

69 

231 

11:60 

3:57 

Location of ®, 

30 Points 

to ® 





5 11 

86 

34 

57 

79 

1:33 

1:19 

10 32 

125 

36 

61 

101 

2:11 

2:01 

20 12 

225 

37 

64 

125 

4:01 

2:04 



TABLE 5 (CONCLUDED) 



Case 

IT 

“cs 

"o 

IT 

“90 

IT 

CIR 

''Tcs 

C^PSI 

Location 

of *, 60 

Points to 

00 





10 

7 

239 

A2 

75 

188 

7:5A 

3:A5 

20 

8 

299 

39 

71 

208 

9:58 

3:5A 

AO 

9 

37A 

37 

69 

231 

11:60 

3:57 

80 

10 

' A63 

36 

68 

257 

15:3A 

A:1A 

Coordinate System 

Convergence Criteria 





17 

52 

85 

97 

lAO 

1:02 

2:0A 

10-3 

18 

88 

71 

92 

lAl 

1:32 

1:55 

10-A 

19 

125 

69 

92 

lAl 

2:11 

1:5A 

10-5 

20 

162 

68 

92 

lAl 

2:A0 

2:11 

Stream Function Convergence 

Criteria 





10-2 

21 

162 

3 

12 

22 

2:A0 

0:35 

10-3 

22 

162 

10 

31 

61 

2:A0 

0:53 

10-A 

23 

162 

35 

61 

101 

2:A0 

1:31 

10-5 

20 

162 

68 

92 

lAl 

2:A0 

2:11 

Note 

lUl: Case numbers 

correspond to 

those of Table 

1. 



Note (?2: CPU time is given in MINUTES rSECOtros. These times are 

subject to some variation depending on the operating 
conditions at the time of the run. 

Legend: 


IT 


CS 


IT, 


IT 


90 


IT 


CIR 
CT, 


CS 


CT 


PS I 


Iterations for coordinate system 

Iterations for 0“ potential flow solution 

Iterations for 90“ potential flow solution 

Iterations for circulation potential flow solution 

CPU time for coordinate system 

CPU time for complete potential flow solution 


TABLE 6 


Coopariion with Analytic Solution - Zero Lift and Other Angles 





Case 




«2 

20.8596 

(zero 

lift) 

33 

-0.0440 

-0.0334 

2.01-01 

8.89-03 

1.6746 

(zero 

lift) 

34 

0.0000 

0.0020 

6.38-03 

5.98-04 

1.65 



35 

0.0000 

-0.0037 

6.12-03 

2.98-03 

2.66 



36 

0.0128 

-0.0358 

1.09-02 

6.87-04 


Legend: 


See Table 2 or 3. Case 33 is Airfoil 111, and the others are Airfoil 

» 2 . 


TABLE 7 


Iterations and Computer Time - Zero Lift and Other Angles 





Case 

IT 

“cs 

'^0 

'^90 

IT.t^ 

'*CS 

^^’psi 

« - -20.8596 

(zero 

lift) 

33 

401 

197 

195 

348 

11:57 

7:36 

0 - - 1.6746 

(zero 

lift) 

34 

419 

131 

170 

350 

12:56 

6:43 

0 ■ - 1.65 



35 

377 

237 

172 

332 

11:10 

7:38 

0 ■ 2.66 



36 

377 

237 

172 

332 

■ 11:10 

7:38 

Legend: See 

Table 

4 or 5. 

Case 

33 Is 

Airfoil 


and the 

others ^ 

are 


Mrfoil » 2 . 


Points OQ 



Tsbls I 



t' • f« . • 

>■» M O « 












TABU 9 (OOMTINUED) 



oriqinal paqe is 

OP POOR QUALITY 




Tablt 9 (COMCLUDBO) OMQfflAL PAGE IS 

09 ! POOR QUALITY. 




TABLE 10 


Karman-Tref ftz Transformation Parameters (Ref. 1 ) 


Transformation: ■= ; — 



k - 1 , 2 , 


\ " 



N 


This transforms a circle in the plane into the airfoil in the plane 
through N successive transformations. 


Airfoil n 
N - 3 

- 0.727277 + iO. 6755902 

- 0 + iO 

- 160“ 

R 2 • 0.39875 - iO. 91706 
« 0 + iO 
^2 • 200 “ 

R 3 “ 0.93667 - iO. 35021 
= -0. 843 - iO. 315866 
“ 12 " 

Circle center = 0 + iO 
Circle radius =1.0 


Airfoil in 
N - 1 

R^ - 1.0 + iO 
- -1.0 + iO , 



Circle Center: -0.04405 + i0.03 

Circle Radius: 1.04448 

Airfoil Chord: 3.9081497 


Airfoil Chord = 3.4647 


TABLE 11 


Pk(ZT 

JF POOR QUALITY 

Changes In Force Coefficients with Increase from Adequate 
Parameter Values ~ Airfoil //I 


Parameter 

Change 

fL 

fa 

Sumber of Points on Airfoil 

37 to 67 

0.009 to 0.002 

0.01 to O.OOA 

dumber of Lines Around Airfoil 

30 to 60 

0,008 to 0.003 

0.004 to 0.003 

Eladius of Outer Boundary 

10 to 20 

0.008 to 0.006 

0.004 to 0.001 

Coordinate Convergence Criteria 

10-*’ to 10-5 

0.005 to 0.002 

no change 

Stream function Convergence Criteria 

10-** to 10-5 

0.01 to 0.002 

0.005 to 0.004 

Kutta Extrapolation Order 

2 to 3 

0.003 to 0.002 

no change 

Kutta Velocity Order 

2 to 3 

worse 

worse 

Velocity Order 

2 to 3 

worse 

0.005 to 0.002 

Pressure Integration Order 

1 to 2 

worse 

worse 


CM 





2000 6 0.7129 

10 0.7642 

20 1.0704 
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(40 lines) 
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Transformed Plane 


Figure 1. 


Field Transformation - Single Body 
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Physical Platie 



Transformed Plane 


Figure 2. Field Transt ormation - Multiple Bodies 


b ) TVo Body 


Figure 3. ConpuCaClonal Grids - Single and Two Body Regions 






Figure 4 


Extrapolacion Points for Application of Kutta Condition 
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Figure 5 . Karman-Tref f tz Airfoil #1 Contours 
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Figure 5. (Concluded) 
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Figure 6. Karman-Tref f tz Airfoil 112 Contours 
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Figure 6. (Concluded) 
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(b.) 67 Points (Case 19) 



(d.) 127 Points (Case 1) 


Figure 8. Effect of Contour Point Distribution on Pressure Distribution - 
Airfoil //I. (Solid line is analytical; symbols are numerical.) 



(a.) 37 Points (Case 16) 
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(b.) 67 Points (Case 19) 


of Contour Point Dititribution on Streamlines - 
line is analytical; cymbols are numerical.) 


Figure 9. Effect 
(Solid 


Airfoil It I 
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(a.) -» dt 5, 30 Points (Case 11) (b.) ® at 10, 60 Points (Case 7) 



(c.) “ at 10, 30 Points (Case 32) (d.) “ at 20, 60 Points (Case 8) 


Figure 11. Effect of Distance to Outer Boundary - Airfoil //I. 

(Solid line is analytical; symbols are numerical.) 
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(a.) Coordinate System Convergence (b.) Coordinate System Convergence 
lO"^ (Case 14) 10‘^ (Case 15) 



(c.) Stream Function Convergence 
10“^ (Case 22) 


(d.) Stream Function Convergence 
10“^ (Case 23) 


Figure 12. Effect of Convergence Criteria - Airfoil ill. (Solid line is 
analytical; symbols are numerical.) 
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(a.) Stagnation Point at Trailing 
Edge (Case 26) 


(b.) First-Order Extrapolation at 
Trailing Edge (Case 24} 




I 


I 

1 

I 

i 

I 

\ 

I 

I 

I 

I 

I 

I 




oiiS" 


•9.96 O.H 




(c.) First-Order Velocity (Case 27) (d.) Third-Order Velocity (Case 28) 


Figure 13. Effect of Extrapolation and Velocity Orders for Kutta 
Condition - Airfoil I/l. (Solid line is analytical; 
symbols are numerical.) 
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(a.) First Order (Case 29) 


l.OOw 


^.? 0 _ 


-s.ooL 


j 

-u.oo' 

•i.og 



I 1 . : : 1 1 

% -g.^o -6.H -o.ug o.»’s o.‘o o.?5 i.so 

1 


(b.) Third Order (Case 30) 

Figure 14. Effect of Order of Velocity in Surface Pressure Calculation - 

Airfoil //I. (Solid line is analytical; symbols are numerical.) 


oF fooh quality 



(a.) 61 Points (Case 6) 


(b.) 67 Points (Case 20) 


•4»mc 

’ ; 


/ 1 


1 

1 I 


i < 

*a9, 3«. 

1 

■f ! 


1 1 - 
4 




1 

* 

-M 31. * 

t 

• 1 


0 



•!.60 •.). 

•d.V? -O.n -O.M v.M 9.W O.n '.»a 

1 
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Figure 

15. Effect of Contour Point 
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(d.) 127 Points (Case 1) 
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Figure 18. Llebeck Laminar Airfoil 








PRESSURE COEFFICIENT 
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Figure 20* Comparison of Experimental and Numerical Pressure 
Distribution - Liebeck Laminar Airfoil 


LIEBECK LAMINAR AIRFOIL 


NUMERICAL 

« TEST DATA[ 22] 



angle of ATTACK 


Figure 21. Comparison of Experimental and Numerical L 
Lift Curve - Liebeck Laminar Airfoil 





(b;) Configuration 



(c.) Configuration 1^2 

Figure 22. Two-Body Coordinate System Configurations. 






























(b.) 61 Points on Each Body 


Figure 28. Coordinate System and Pressure Distribution for 
Two Circles, Configuration 112 . 


















